﻿#region Copyright information
// 
// Copyright © 2005-2013 Yongkee Cho. All rights reserved.
// 
// This code is a part of the Biological Object Library and governed under the terms of the
// GNU Lesser General  Public License (LGPL) version 2.1 which accompanies this distribution.
// For more information on the LGPL, please visit http://bol.codeplex.com/license.
// 
// - Filename: WeibullDistribution.cs
// - Author: Yongkee Cho
// - Email: yongkeecho@gmail.com
// - Date Created: 2012-09-06 11:39 AM
// - Last Modified: 2013-01-25 3:59 PM
// 
#endregion
using System;
using System.Collections.Generic;
using BOL.Maths.Functions;

namespace BOL.Maths.Distributions
{
    /// <summary>
    /// Represents a Weibull with lambda scale factor and k shape factor.
    /// <see>Mathematica: ref/WeibullDistribution</see>
    /// <see>http://en.wikipedia.org/wiki/Weibull_distribution</see>
    /// </summary>
    public sealed class WeibullDistribution : IUnivariateDistribution<double>, IEquatable<WeibullDistribution>
    {
        #region Private variables

        private double _lambda; // scale factor
        private double _k; // shape factor
        private double _oneOverK;
        private const double Gamma = 0.57721566490153286060651209008240243104215933593992;

        #endregion

        #region Public properties

        /// <summary>Gets or sets the scale factor of a Weibull distribution.</summary>
        public double Lambda { get { return _lambda; } set { if (value <= 0) throw new ArgumentOutOfRangeException("value", "lambda is expected to be a positive real number."); _lambda = value; } }

        /// <summary>Gets or sets the shape factor of freedom of a Weibull distribution.</summary>
        public double K { get { return _k; } set { if (value <= 0) throw new ArgumentOutOfRangeException("value", "k is expected to be a positive real number."); _k = value; _oneOverK = 1 / value; } }

        /// <summary>Gets the range of a Weibull distribution.</summary>
        public IRange<double> Domain { get { return new Range<double>(0, double.PositiveInfinity); } }

        /// <summary>Gets the mean of a Weibull distribution.</summary>
        public double Mean { get { return _lambda * DistributionFunctions.Gamma(1 + _oneOverK); } }

        /// <summary>Gets the median of a Weibull distribution.</summary>
        public double Median { get { return _lambda * Math.Pow(Math.Log(2), _oneOverK); } }

        /// <summary>Gets the mode of a Weibull distribution.</summary>
        public double Mode { get { return _lambda * Math.Pow((_k - 1) / _k, _oneOverK); } }

        /// <summary>Gets the variance of a Weibull distribution.</summary>
        public double Variance { get { return _lambda * _lambda * DistributionFunctions.Gamma(1 + 2 * _oneOverK) - Math.Pow(Mean, 2); } }

        /// <summary>Gets the skewness of a Weibull distribution.</summary>
        public double Skewness { get { return (DistributionFunctions.Gamma(1 + 3 / _k) * _lambda * _lambda * _lambda - 3 * Mean * Variance - Math.Pow(Mean, 3)) / Math.Pow(Variance, 3 / 2); } }

        /// <summary>Gets the kurtosis of a Weibull distribution.</summary>
        public double Kurtosis { get { throw new NotImplementedException(); } }

        /// <summary>Gets the entropy of a Weibull distribution.</summary>
        public double Entropy { get { return Gamma * (1 - _oneOverK) + Math.Log(_lambda * _oneOverK) + 1; } }

        #endregion

        #region Constructor

        /// <summary>
        /// Instantiates a Weibull distribution.
        /// </summary>
        /// <param name="lambda">scale</param>
        /// <param name="k">shape</param>
        public WeibullDistribution(double lambda, double k)
        {
            if (lambda <= 0)
                throw new ArgumentOutOfRangeException("lambda", "lambda is expected to be a positive real number.", "lambda");
            if (_k <= 0)
                throw new ArgumentOutOfRangeException("k", "k is expected to be a positive real number.", "k");

            _lambda = lambda;
            _k = k;
            _oneOverK = 1 / _k;
        }

        #endregion

        #region ICloneable implementation

        public WeibullDistribution Clone()
        {
            return new WeibullDistribution(_lambda, _k);
        }

        object ICloneable.Clone()
        {
            return Clone();
        }
        
        #endregion

        #region Public methods

        /// <summary>
        /// Returns the Weibull density given value.
        /// </summary>
        /// <param name="value"></param>
        /// <returns>Weibull density</returns>
        public double Pdf(double value)
        {
            var temp = value / _lambda;
            return value >= 0 ? (_k / _lambda) * Math.Pow(temp, _k - 1) * Math.Exp(-Math.Pow(temp, _k)) : 0;
        }

        public double Hazard(double value)
        {
            return value >= 0 ? (_k / _lambda) * Math.Pow(value / _lambda, _k - 1) : 0;
        }

        /// <summary>
        /// Returns the cummulative Weibull density given value.
        /// </summary>
        /// <param name="value"></param>
        /// <returns>cummulative Weibull density</returns>
        public double Cdf(double value)
        {
            return 1 - Math.Exp(- Math.Pow(value / _lambda, _k));
        }

        public double CummulativeHazard(double value)
        {
            return value >= 0 ? Math.Pow(value / _lambda, _k) : 0;
        }

        public double Survival(double value)
        {
            return value >= 0 ? Math.E - Math.Pow(value / _lambda, _k) : 0;
        }

        public double InverseSurvival(double p)
        {
            return (p.Equals(1.0) || p < 0) ? 0 : _lambda * Math.Pow(-Math.Log(p), 1 / _k);
        }

        /// <summary>
        /// Returns the value of the Weibull distribution for given probability.
        /// </summary>
        /// <param name="p">probability</param>
        /// <returns></returns>
        public double Quantile(double p)
        {
            return (p.Equals(1.0) || p < 0) ? 0 : _lambda * Math.Pow(-Math.Log(1 - p), 1 / _k);
        }

        /// <summary>
        /// Returns a random value from the Weibull distribution.
        /// </summary>
        /// <returns>random value</returns>
        public double Sample(Random r)
        {
            throw new NotImplementedException();
        }

        /// <summary>
        /// Estimates parameters of the Weibull distribution using maximum likelihood.
        /// </summary>
        public void MaximumLikelihoodEstimate(IEnumerable<double> source)
        {
            throw new NotImplementedException();
        }

        /// <summary>
        /// Estimates parameters of the Weibull distribution using Bayesian methods.
        /// </summary>
        public void BayesianEstimate(IEnumerable<double> source, params dynamic[] priorParameters)
        {
            throw new NotImplementedException();
        }

        #endregion

        #region IEquatable<T> implementation

        public bool Equals(WeibullDistribution other)
        {
            return _lambda.Equals(other._lambda) & _k.Equals(other._k);
        }

        #endregion

        #region Object overriden

        public override int GetHashCode()
        {
            return Lambda.GetHashCode() ^ K.GetHashCode();
        }

        public override bool Equals(object other)
        {
            if (other == null)
                throw new ArgumentNullException("other");

            if (!(other is WeibullDistribution))
                throw new InvalidCastException("The 'other' argument is not a WeibullDistribution object.");
            
            return Equals(other as WeibullDistribution);
        }

        public override string ToString()
        {
            return String.Format("Weibull(lambda = {0:G4}, k = {1:G4})", _lambda, _k);
        }

        #endregion
    }
}
